DISTRIBUTIONAL EQUITY

Reference categories: We will chose those with more responses (migrant1, females, married, secondary education). We only have to change the reference category for Education, because the other ones are already set.

data_fiji_dist<- within(data_fiji_dist, Education_cat3 <- relevel(Education_cat3, ref = 2)) #set secondary education as the reference category as it has more counts

data_fiji_dist$Distr_Eq_int<-data_fiji_dist$Distr_Eq_Impact_int

data_fiji_dist$Distr_Eq_int_cat<-as.factor(data_fiji_dist$Distr_Eq_int)

Bayesian (brms package)

The cumulative model assumes that the effect of a predictor is equal across the rating categories (e.g. the effect of marital status is the same across the levels of distributional equity). However, this assumption may not be appropriate, and marital status may impact rating categories differently. If the effects of predictors vary across rating categories, we call the resulting model to have category specific (cs) effects.

Checking proportional odds assumption with brms

(https://discourse.mc-stan.org/t/proportional-odds-assumptions/5316)

We need to investigate if each predictor has category specific effects. In order to do this we use an adjacent category model which uses the family=acat() instead of family=cumulative() (Burkner and Vuorre xx Ordinal regression models in psychology: a tutorial https://psyarxiv.com/x8swp/).

Additive model

#M.de5_add<-brm(Distr_Eq_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ #(1|village),data=data_fiji_dist,family=cumulative("logit"),control = list(adapt_delta = 0.99))  

#saveRDS(M.de5_add,file="M.de5_add.rda")

M.de5_add<-readRDS(file="M.de5_add.rda")

prior_summary(M.de5_add)# check default priors

##                 prior     class            coef   group resp dpar nlpar lb ub
##                (flat)         b                                              
##                (flat)         b      Age_num_sd                              
##                (flat)         b Education_cat31                              
##                (flat)         b Education_cat33                              
##                (flat)         b         Genderm                              
##                (flat)         b Marital_statuss                              
##                (flat)         b Marital_statusw                              
##                (flat)         b       migrant22                              
##                (flat)         b      wealth2_sd                              
##  student_t(3, 0, 2.5) Intercept                                              
##  student_t(3, 0, 2.5) Intercept               1                              
##  student_t(3, 0, 2.5) Intercept               2                              
##  student_t(3, 0, 2.5) Intercept               3                              
##  student_t(3, 0, 2.5) Intercept               4                              
##  student_t(3, 0, 2.5)        sd                                          0   
##  student_t(3, 0, 2.5)        sd                 village                  0   
##  student_t(3, 0, 2.5)        sd       Intercept village                  0   
##        source
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##  (vectorized)

Add weakly informative priors.

Weakly informative priors

priors2 <-prior(normal(0, 2.5), class = "b")+
  prior(normal(0, 2.5), class = "Intercept")+
  prior(student_t(3,0,2), class = "sd")

Modelling cs(Education_cat3) does not work.We will model it without specific categories for education

#```{r } M.de5_acat_all_noed<-brm(bf(Distr_Eq_int~cs(Gender)+cs(migrant2)+Education_cat3+cs(wealth2_sd)+cs(Age_num_sd)+cs(Marital_status)+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=acat(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.de5_acat_all_noed,file="M.de5_acat_all_noed.rda")

M.de5_acat_all_noed<-readRDS(file="M.de5_acat_all_noed.rda")

loo.M.de5_acat_all_noed<-loo(M.de5_acat_all_noed)

## Warning: Found 2 observations with a pareto_k > 0.7 in model
## 'M.de5_acat_all_noed'. It is recommended to set 'moment_match = TRUE' in order
## to perform moment matching for problematic observations.

M.de5_acat_all2<-readRDS(file="M.de5_acat_all2.rda")

loo.M.de5_acat_all2<-loo(M.de5_acat_all2)

#```{r } M.de5_acat_all2<-brm(bf(Distr_Eq_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=acat(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.de5_acat_all2,file="M.de5_acat_all2.rda")

M.de5_acat_all2<-readRDS(file="M.de5_acat_all2.rda")

loo.M.de5_acat_all2<-loo(M.de5_acat_all2)

Additive model

#```{r additive model2, cashe=TRUE} M.de5_add_2<-brm(bf(Distr_Eq_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#{r} saveRDS(M.de5_add_2,file="M.de5_add_2.rda") #

M.de5_add_2<-readRDS(file="M.de5_add_2.rda")

pars <- M.de5_add_2%>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.de5_add_2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.de5_add_2)

 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)

 y_int <- data_fiji_dist$Distr_Eq_int

ppc_bars(y_int, yrep_int)

Trace plots:

M.de5_add_2%>%mcmc_trace

Chain convergence:

M.de5_add_2%>% mcmc_plot(type='rhat_hist')

## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_add_2%>% mcmc_plot(type='neff_hist')

## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.de5_add_2%>% mcmc_plot(type='acf_bar')

## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.de5_add_2<- as.array(M.de5_add_2)

Residual plot:

preds<-posterior_predict(M.de5_add_2,ndraws=250,summary=FALSE)

M.de5_add_2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_dist$Distr_Eq_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.de5_add_2.resids) #residuals look good

summary(M.de5_add_2)

##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Distr_Eq_int ~ Gender + migrant2 + Education_cat3 + wealth2_sd + Age_num_sd + Marital_status + (1 | village) 
##    Data: data_fiji_dist (Number of observations: 159) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.26     0.02     1.02 1.00     2623     3148
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]       -4.31      0.58    -5.48    -3.24 1.00     3213     3217
## Intercept[2]       -2.40      0.44    -3.31    -1.55 1.00     3150     3174
## Intercept[3]       -1.82      0.43    -2.69    -1.02 1.00     3209     3027
## Intercept[4]        0.57      0.40    -0.22     1.37 1.00     3105     3089
## Genderm             0.09      0.37    -0.64     0.82 1.00     2960     3180
## migrant22          -0.71      0.35    -1.39    -0.03 1.00     3277     3008
## Education_cat31    -0.85      0.42    -1.66    -0.04 1.00     3208     3173
## Education_cat33    -0.99      0.40    -1.78    -0.19 1.00     3135     3093
## wealth2_sd          0.17      0.34    -0.52     0.82 1.00     3069     2976
## Age_num_sd          0.50      0.39    -0.26     1.29 1.00     3307     3174
## Marital_statuss     0.09      0.52    -0.97     1.11 1.00     3246     3172
## Marital_statusw    -1.33      0.52    -2.38    -0.32 1.00     3052     3071
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

loo.M.de5_add_2<-loo(M.de5_add_2)

loo.M.de5_add_2

## 
## Computed from 3200 by 159 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -212.4  9.8
## p_loo        16.3  1.3
## looic       424.9 19.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Compare cumulative and acat

loo_compare(loo.M.de5_add_2,loo.M.de5_acat_all2,loo.M.de5_acat_all_noed)

##                     elpd_diff se_diff
## M.de5_add_2           0.0       0.0  
## M.de5_acat_all2      -1.3       1.5  
## M.de5_acat_all_noed -11.6       4.9

Gender x Migrant

#```{r genderxmigrant model2, cashe=TRUE} M.de5_int.mig2<-brm(bf(Distr_Eq_int~Gender*migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#{r} saveRDS(M.de5_int.mig2,file="M.de5_int.mig2.rda") #

M.de5_int.mig2<-readRDS(file="M.de5_int.mig2.rda")

pars <- M.de5_int.mig2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.de5_int.mig2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.de5_int.mig2)

 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)

 y_int <- data_fiji_dist$Distr_Eq_int

ppc_bars(y_int, yrep_int)

Trace plots:

M.de5_int.mig2%>%mcmc_trace

Chain convergence:

M.de5_int.mig2 %>% mcmc_plot(type='rhat_hist')

## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_int.mig2 %>% mcmc_plot(type='neff_hist')

## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.de5_int.mig2 %>% mcmc_plot(type='acf_bar')

## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.de5_int.mig2<- as.array(M.de5_int.mig2)

mcmc_acf(posterior_M.de5_int.mig2, 
         pars = c("b_Intercept[1]","b_Intercept[2]","b_Intercept[3]","b_Intercept[4]", "b_Genderm", "b_migrant22", "b_Education_cat31","b_Education_cat33", "b_wealth2_sd","b_Age_num_sd","b_Marital_statuss","b_Marital_statusw"))

Residual plot:

preds<-posterior_predict(M.de5_int.mig2,ndraws=250,summary=FALSE)

M.de5_int.mig2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_dist$Distr_Eq_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.de5_int.mig2.resids) #residuals look good

summary(M.de5_int.mig2)

##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Distr_Eq_int ~ Gender * migrant2 + Education_cat3 + wealth2_sd + Age_num_sd + Marital_status + (1 | village) 
##    Data: data_fiji_dist (Number of observations: 159) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.43      0.26     0.03     1.04 1.00     2558     2909
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]         -4.25      0.59    -5.43    -3.16 1.00     3084     3052
## Intercept[2]         -2.30      0.45    -3.22    -1.44 1.00     3215     3197
## Intercept[3]         -1.70      0.43    -2.58    -0.87 1.00     3205     2963
## Intercept[4]          0.73      0.42    -0.10     1.56 1.00     3339     3116
## Genderm               0.39      0.42    -0.45     1.22 1.00     3209     3131
## migrant22            -0.27      0.43    -1.13     0.59 1.00     2985     3161
## Education_cat31      -0.90      0.42    -1.71    -0.07 1.00     3102     3175
## Education_cat33      -1.06      0.39    -1.86    -0.32 1.00     2924     3183
## wealth2_sd            0.20      0.33    -0.45     0.86 1.00     3138     3091
## Age_num_sd            0.59      0.40    -0.19     1.37 1.00     3188     3120
## Marital_statuss       0.13      0.52    -0.88     1.16 1.00     3275     3153
## Marital_statusw      -1.34      0.54    -2.40    -0.29 1.00     3312     3127
## Genderm:migrant22    -1.19      0.70    -2.55     0.18 1.00     3112     3299
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

mcmc_intervals(posterior_M.de5_int.mig2, prob = 0.8, # 80% intervals
  prob_outer = 0.95, # 95%
  point_est = "mean"
,pars = c("b_Genderm", "b_migrant22", "b_Education_cat31","b_Education_cat33", "b_wealth2_sd","b_Age_num_sd","b_Marital_statuss","b_Marital_statusw","b_Genderm:migrant22"))

Conditional effects:

Ploting the interaction

conditions <- data.frame(Gender = c("f","m"))

conditional_effects(M.de5_int.mig2, "migrant2", prob=0.95,
                        categorical=TRUE, conditions=conditions) #re_formula=NULL takes into account random effects

inter_plot<-conditional_effects(M.de5_int.mig2, "migrant2", prob=0.95,
                        categorical=TRUE, conditions=conditions) #re_formula=NULL takes into account random effects

plot1<- plot(inter_plot, plot = FALSE)[[1]] + ##ggplot object
  theme_classic()+scale_color_manual(values = c("1"= "darkred", 
                             "2"= "red3", 
                             "3"= "grey61",
                             "4" = "steelblue2",
                             "5" = "royalblue4")) +
  labs(x="Migrant Status")+
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=15),
        legend.text=element_text(size=10),
        legend.title=element_text(size=10))+
  scale_fill_discrete(labels = c("Very unfair", "Unfair", "Neutral","Fair", "Very fair"))

plot1

inter_plot.80<-conditional_effects(M.de5_int.mig2, "migrant2", prob=0.80,
                        categorical=TRUE, conditions=conditions) #re_formula=NULL takes into account random effects

plot2<- plot(inter_plot.80, plot = FALSE)[[1]] + ##ggplot object
  theme_classic()+scale_color_manual(values = c("1"= "darkred", 
                             "2"= "red3", 
                             "3"= "grey61",
                             "4" = "steelblue2",
                             "5" = "royalblue4")) +
  labs(x="Migrant Status")+
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=15),
        legend.text=element_text(size=10),
        legend.title=element_text(size=10))

plot2

Change colors

plot1

plot1<- plot(inter_plot, plot = FALSE)[[1]] + ##ggplot object
  theme_classic()+scale_color_manual(values = c("1"= "#CC0033", 
                             "2"= "#FF6699", 
                             "3"= "#CCCCCC",
                             "4" = "#CCCCFF",
                             "5" = "#6699FF")) +
  labs(x="Migrant Status")+
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=15),
        legend.text=element_text(size=10),
        legend.title=element_text(size=10))

plot1

loo.M.de5_int.mig2<-loo(M.de5_int.mig2)

loo.M.de5_int.mig2

## 
## Computed from 3200 by 159 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -211.7 10.1
## p_loo        17.4  1.5
## looic       423.4 20.2
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

M.de5_int.marital3<-readRDS(file="M.de5_int.marital3.rda")

M.de5_int.marital2<-M.de5_int.marital3

pars <- M.de5_int.marital2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.de5_int.marital2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

 yrep_char <- posterior_predict(M.de5_int.marital2)

 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)

 y_int <- data_fiji_dist$Distr_Eq_int

ppc_bars(y_int, yrep_int)

M.de5_int.marital2%>%mcmc_trace

M.de5_int.marital2%>% mcmc_plot(type='rhat_hist')

## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_int.marital2%>% mcmc_plot(type='neff_hist')

## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_int.marital2%>% mcmc_plot(type='acf_bar')

## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.de5_int.marital2<- as.array(M.de5_int.marital2)

preds<-posterior_predict(M.de5_int.marital2,ndraws=250,summary=FALSE)

M.de5_int.marital2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_dist$Distr_Eq_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.de5_int.marital2.resids) #residuals look good

summary(M.de5_int.marital2)

##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Distr_Eq_int ~ Gender * Marital_status + Education_cat3 + wealth2_sd + Age_num_sd + migrant2 + (1 | village) 
##    Data: data_fiji_dist (Number of observations: 159) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.41      0.26     0.03     1.01 1.00     2371     2262
## 
## Population-Level Effects: 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]               -4.29      0.62    -5.57    -3.15 1.00     3287
## Intercept[2]               -2.37      0.47    -3.32    -1.47 1.00     3093
## Intercept[3]               -1.79      0.45    -2.70    -0.91 1.00     3180
## Intercept[4]                0.60      0.43    -0.21     1.46 1.00     3074
## Genderm                     0.10      0.42    -0.70     0.92 1.00     3175
## Marital_statuss             0.10      0.71    -1.30     1.47 1.00     2994
## Marital_statusw            -1.34      0.60    -2.53    -0.18 1.00     3119
## Education_cat31            -0.74      0.40    -1.54     0.03 1.00     3178
## Education_cat33            -0.97      0.39    -1.73    -0.20 1.00     3128
## wealth2_sd                  0.16      0.34    -0.51     0.83 1.00     2859
## Age_num_sd                  0.46      0.39    -0.29     1.23 1.00     3055
## migrant22                  -0.70      0.36    -1.39     0.00 1.00     3015
## Genderm:Marital_statuss    -0.05      0.95    -1.93     1.73 1.00     2974
## Genderm:Marital_statusw     0.08      0.92    -1.69     1.86 1.00     2526
##                         Tail_ESS
## Intercept[1]                2831
## Intercept[2]                3057
## Intercept[3]                3137
## Intercept[4]                2962
## Genderm                     3053
## Marital_statuss             3067
## Marital_statusw             2859
## Education_cat31             3023
## Education_cat33             3115
## wealth2_sd                  3027
## Age_num_sd                  3039
## migrant22                   2907
## Genderm:Marital_statuss     3158
## Genderm:Marital_statusw     3081
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

loo.M.de5_int.marital2<-loo(M.de5_int.marital2)

loo.M.de5_int.marital2

## 
## Computed from 3200 by 159 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -215.0  9.8
## p_loo        18.4  1.6
## looic       430.1 19.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x education

#```{r M.de5_int.edu2, cashe=TRUE} M.de5_int.edu2<-brm(bf(Distr_Eq_int~Gender*Education_cat3+Marital_status+wealth2_sd+Age_num_sd+migrant2+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.de5_int.edu2,file="M.de5_int.edu2.rda")

M.de5_int.edu2<-readRDS(file="M.de5_int.edu2.rda")

pars <- M.de5_int.edu2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.de5_int.edu2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.de5_int.edu2)

 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)

 y_int <- data_fiji_dist$Distr_Eq_int

ppc_bars(y_int, yrep_int)

Trace plots:

M.de5_int.edu2%>%mcmc_trace

Chain convergence:

M.de5_int.edu2%>% mcmc_plot(type='rhat_hist')

## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_int.edu2%>% mcmc_plot(type='neff_hist')

## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.de5_int.edu2%>% mcmc_plot(type='acf_bar')

## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.de5_int.edu2<- as.array(M.de5_int.edu2)

Residual plot:

preds<-posterior_predict(M.de5_int.edu2,ndraws=250,summary=FALSE)

M.de5_int.edu2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_dist$Distr_Eq_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.de5_int.edu2.resids) #residuals look good

summary(M.de5_int.edu2)

##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Distr_Eq_int ~ Gender * Education_cat3 + Marital_status + wealth2_sd + Age_num_sd + migrant2 + (1 | village) 
##    Data: data_fiji_dist (Number of observations: 159) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.27     0.02     0.98 1.00     2691     2685
## 
## Population-Level Effects: 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]               -4.28      0.59    -5.46    -3.17 1.00     3143
## Intercept[2]               -2.37      0.45    -3.26    -1.50 1.00     2702
## Intercept[3]               -1.79      0.43    -2.64    -0.95 1.00     2868
## Intercept[4]                0.60      0.42    -0.20     1.43 1.00     2926
## Genderm                     0.09      0.45    -0.77     0.99 1.00     3296
## Education_cat31            -0.84      0.58    -1.98     0.33 1.00     3286
## Education_cat33            -0.87      0.53    -1.89     0.13 1.00     3036
## Marital_statuss             0.09      0.54    -0.94     1.15 1.00     3022
## Marital_statusw            -1.30      0.54    -2.37    -0.26 1.00     3154
## wealth2_sd                  0.18      0.34    -0.51     0.87 1.00     3067
## Age_num_sd                  0.48      0.40    -0.31     1.26 1.00     3496
## migrant22                  -0.72      0.35    -1.40    -0.03 1.00     3202
## Genderm:Education_cat31     0.20      0.77    -1.31     1.72 1.00     3180
## Genderm:Education_cat33    -0.18      0.74    -1.58     1.24 1.00     3042
##                         Tail_ESS
## Intercept[1]                3303
## Intercept[2]                2992
## Intercept[3]                2988
## Intercept[4]                2848
## Genderm                     3217
## Education_cat31             3221
## Education_cat33             3258
## Marital_statuss             2997
## Marital_statusw             3194
## wealth2_sd                  2969
## Age_num_sd                  2895
## migrant22                   3174
## Genderm:Education_cat31     3112
## Genderm:Education_cat33     3010
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

loo.M.de5_int.edu2<-loo(M.de5_int.edu2)

loo.M.de5_int.edu2

## 
## Computed from 3200 by 159 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -215.1  9.9
## p_loo        18.5  1.5
## looic       430.2 19.7
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x age

#```{r M.de5_int.age2, cashe=TRUE} M.de5_int.age2<-brm(bf(Distr_Eq_int~Gender*Age_num_sd+Education_cat3+Marital_status+wealth2_sd+migrant2+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#{r} saveRDS(M.de5_int.age2,file="M.de5_int.age2.rda") #

M.de5_int.age2<-readRDS(file="M.de5_int.age2.rda")

pars <- M.de5_int.age2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.de5_int.age2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.de5_int.age2)

 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)

 y_int <- data_fiji_dist$Distr_Eq_int

ppc_bars(y_int, yrep_int)

Trace plots:

M.de5_int.age2%>%mcmc_trace

Chain convergence:

M.de5_int.age2%>% mcmc_plot(type='rhat_hist')

## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_int.age2%>% mcmc_plot(type='neff_hist')

## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.de5_int.age2%>% mcmc_plot(type='acf_bar')

## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.de5_int.age2<- as.array(M.de5_int.age2)

Residual plot:

preds<-posterior_predict(M.de5_int.age2,ndraws=250,summary=FALSE)

M.de5_int.age2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_dist$Distr_Eq_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.de5_int.age2.resids) #residuals look good

summary(M.de5_int.age2)

##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Distr_Eq_int ~ Gender * Age_num_sd + Education_cat3 + Marital_status + wealth2_sd + migrant2 + (1 | village) 
##    Data: data_fiji_dist (Number of observations: 159) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.41      0.26     0.03     1.03 1.00     2807     3058
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]          -4.31      0.58    -5.46    -3.22 1.00     2966     3173
## Intercept[2]          -2.40      0.44    -3.27    -1.54 1.00     3264     3144
## Intercept[3]          -1.81      0.43    -2.67    -0.99 1.00     3191     3061
## Intercept[4]           0.58      0.41    -0.23     1.40 1.00     3319     2529
## Genderm                0.09      0.37    -0.60     0.84 1.00     3289     3127
## Age_num_sd             0.45      0.57    -0.64     1.55 1.00     2917     2972
## Education_cat31       -0.85      0.42    -1.67    -0.03 1.00     2899     3215
## Education_cat33       -0.99      0.40    -1.78    -0.21 1.00     3189     3116
## Marital_statuss        0.08      0.52    -0.96     1.07 1.00     2873     3089
## Marital_statusw       -1.31      0.55    -2.39    -0.23 1.00     3052     2710
## wealth2_sd             0.17      0.35    -0.51     0.85 1.00     3329     2982
## migrant22             -0.72      0.36    -1.44    -0.04 1.00     3363     3172
## Genderm:Age_num_sd     0.08      0.66    -1.21     1.39 1.00     2693     3227
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

loo.M.de5_int.age2<-loo(M.de5_int.age2)

loo.M.de5_int.age2

## 
## Computed from 3200 by 159 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -213.4  9.9
## p_loo        17.2  1.4
## looic       426.7 19.8
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x wealth

#```{r M.de5_int.wealth2, cashe=TRUE} M.de5_int.wealth2<-brm(bf(Distr_Eq_int~Gender*wealth2_sd+Age_num_sd+Education_cat3+Marital_status+migrant2+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.de5_int.wealth2,file="M.de5_int.wealth2.rda")

M.de5_int.wealth2<-readRDS(file="M.de5_int.wealth2.rda")

pars <- M.de5_int.wealth2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.de5_int.wealth2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.de5_int.wealth2)

 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)

 y_int <- data_fiji_dist$Distr_Eq_int

ppc_bars(y_int, yrep_int)

Trace plots:

M.de5_int.wealth2%>%mcmc_trace

Chain convergence:

M.de5_int.wealth2%>% mcmc_plot(type='rhat_hist')

## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_int.wealth2%>% mcmc_plot(type='neff_hist')

## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.de5_int.wealth2%>% mcmc_plot(type='acf_bar')

## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.de5_int.wealth2<- as.array(M.de5_int.wealth2)

Residual plot:

preds<-posterior_predict(M.de5_int.wealth2,ndraws=250,summary=FALSE)

M.de5_int.wealth2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_dist$Distr_Eq_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.de5_int.wealth2.resids) #residuals look good

summary(M.de5_int.wealth2)

##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Distr_Eq_int ~ Gender * wealth2_sd + Age_num_sd + Education_cat3 + Marital_status + migrant2 + (1 | village) 
##    Data: data_fiji_dist (Number of observations: 159) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.37      0.25     0.02     0.97 1.00     2687     2897
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]          -4.35      0.59    -5.61    -3.24 1.00     3063     2965
## Intercept[2]          -2.42      0.44    -3.30    -1.58 1.00     2962     3015
## Intercept[3]          -1.84      0.42    -2.70    -1.02 1.00     2902     2824
## Intercept[4]           0.56      0.40    -0.21     1.33 1.00     2805     3181
## Genderm                0.08      0.36    -0.62     0.77 1.00     3057     3107
## wealth2_sd             0.41      0.41    -0.38     1.21 1.00     2749     2723
## Age_num_sd             0.51      0.39    -0.25     1.27 1.00     2864     3173
## Education_cat31       -0.84      0.41    -1.66    -0.06 1.00     3220     2848
## Education_cat33       -0.97      0.39    -1.76    -0.23 1.00     2974     3117
## Marital_statuss        0.06      0.53    -0.98     1.12 1.00     3056     3303
## Marital_statusw       -1.37      0.53    -2.39    -0.35 1.00     2755     2987
## migrant22             -0.71      0.35    -1.39    -0.04 1.00     3083     3134
## Genderm:wealth2_sd    -0.63      0.65    -1.91     0.68 1.00     3139     2822
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

loo.M.de5_int.wealth2<-loo(M.de5_int.wealth2)

loo.M.de5_int.wealth2

## 
## Computed from 3200 by 159 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -212.6  9.9
## p_loo        16.6  1.4
## looic       425.3 19.7
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Compare models

loo_compare(loo.M.de5_int.wealth2,loo.M.de5_int.age2,loo.M.de5_int.edu2,loo.M.de5_int.marital2,loo.M.de5_int.mig2,loo.M.de5_add_2 )

##                    elpd_diff se_diff
## M.de5_int.mig2      0.0       0.0   
## M.de5_add_2        -0.8       2.0   
## M.de5_int.wealth2  -0.9       2.1   
## M.de5_int.age2     -1.7       2.0   
## M.de5_int.marital2 -3.3       2.1   
## M.de5_int.edu2     -3.4       2.1

PROCEDURAL EQUITY

data_fiji_proced<- within(data_fiji_proced, Education_cat3 <- relevel(Education_cat3, ref = 2)) #set secondary education as the reference category as it has more counts

str(data_fiji_proced$Proced_Equity_int)

##  int [1:168] 4 5 5 5 5 4 4 4 4 3 ...

#```{r } M.pe5_acat_all_noed<-brm(bf(Proced_Equity_int~cs(Gender)+cs(migrant2)+cs(Education_cat3)+cs(wealth2_sd)+cs(Age_num_sd)+cs(Marital_status)+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=acat(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.pe5_acat_all_noed,file="M.pe5_acat_all_noed.rda")

M.pe5_acat_all_noed<-readRDS(file="M.pe5_acat_all_noed.rda")

M.pe5_acat_all_noed<-loo(M.pe5_acat_all_noed)

#```{r } M.pe5_acat_all2<-brm(bf(Proced_Equity_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=acat(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.pe5_acat_all2,file="M.pe5_acat_all2.rda")

M.pe5_acat_all2<-readRDS(file="M.pe5_acat_all2.rda")

M.pe5_acat_all2<-readRDS(file="M.pe5_acat_all2.rda")

loo.M.pe5_acat_all2<-loo(M.pe5_acat_all2)

Additive model

#```{r proced additive model2, cashe=TRUE} M.pe5_add_2<-brm(bf(Proced_Equity_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.pe5_add_2,file="M.pe5_add_2.rda")

M.pe5_add_2<-readRDS(file="M.pe5_add_2.rda")

pars <- M.pe5_add_2%>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.pe5_add_2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.pe5_add_2)

 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)

 y_int <- data_fiji_proced$Proced_Equity_int

ppc_bars(y_int, yrep_int)

Trace plots:

M.pe5_add_2%>%mcmc_trace

Chain convergence:

M.pe5_add_2%>% mcmc_plot(type='rhat_hist')

## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.pe5_add_2%>% mcmc_plot(type='neff_hist')

## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.pe5_add_2%>% mcmc_plot(type='acf_bar')

## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.pe5_add_2<- as.array(M.pe5_add_2)

Residual plot:

preds<-posterior_predict(M.pe5_add_2,ndraws=250,summary=FALSE)

M.pe5_add_2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_proced$Proced_Equity_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.pe5_add_2.resids) #residuals look good

summary(M.pe5_add_2)

##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Proced_Equity_int ~ Gender + migrant2 + Education_cat3 + wealth2_sd + Age_num_sd + Marital_status + (1 | village) 
##    Data: data_fiji_proced (Number of observations: 168) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.42      0.27     0.02     1.02 1.00     2559     2732
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]       -3.03      0.47    -3.99    -2.16 1.00     3324     3252
## Intercept[2]       -1.86      0.41    -2.70    -1.07 1.00     3356     3172
## Intercept[3]       -1.31      0.40    -2.11    -0.53 1.00     3304     3188
## Intercept[4]        0.63      0.39    -0.11     1.40 1.00     3351     3258
## Genderm             0.31      0.36    -0.40     1.00 1.00     3175     2941
## migrant22          -0.52      0.32    -1.14     0.11 1.00     3435     3256
## Education_cat31    -0.24      0.40    -1.03     0.55 1.00     3234     3161
## Education_cat33    -0.91      0.37    -1.62    -0.19 1.00     3234     3128
## wealth2_sd          0.10      0.31    -0.52     0.71 1.00     2907     2595
## Age_num_sd          0.21      0.35    -0.48     0.89 1.00     3331     3130
## Marital_statuss     0.55      0.50    -0.41     1.52 1.00     3228     3070
## Marital_statusw    -0.23      0.50    -1.15     0.76 1.00     3343     3064
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

loo.M.pe5_add_2<-loo(M.pe5_add_2)

loo.M.pe5_add_2

## 
## Computed from 3200 by 168 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -239.0  9.2
## p_loo        16.1  1.1
## looic       478.1 18.5
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Compare cumulative and acat

loo_compare(loo.M.pe5_add_2,loo.M.pe5_acat_all2,M.pe5_acat_all_noed)

##                     elpd_diff se_diff
## M.pe5_add_2           0.0       0.0  
## M.pe5_acat_all2      -0.9       1.6  
## M.pe5_acat_all_noed -11.8       5.6

posterior_M.pe5_add_2<- as.array(M.pe5_add_2)

mcmc_intervals(posterior_M.pe5_add_2, prob = 0.8, # 80% intervals
  prob_outer = 0.95, # 95%
  point_est = "mean"
,pars = c("b_Genderm", "b_migrant22", "b_Education_cat31","b_Education_cat33", "b_wealth2_sd","b_Age_num_sd","b_Marital_statuss","b_Marital_statusw"))

Gender x Migrant

#```{r proced genderxmigrant model2, cashe=TRUE} M.pe5_int.mig2<-brm(bf(Proced_Equity_int~Gender*migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.pe5_int.mig2,file="M.pe5_int.mig2.rda")

M.pe5_int.mig2<-readRDS(file="M.pe5_int.mig2.rda")

pars <- M.pe5_int.mig2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.pe5_int.mig2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.pe5_int.mig2)

 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)

 y_int <- data_fiji_proced$Proced_Equity_int

ppc_bars(y_int, yrep_int)

Trace plots:

M.pe5_int.mig2%>%mcmc_trace

Chain convergence:

M.pe5_int.mig2 %>% mcmc_plot(type='rhat_hist')

## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.pe5_int.mig2 %>% mcmc_plot(type='neff_hist')

## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.pe5_int.mig2 %>% mcmc_plot(type='acf_bar')

## Warning: Removed 84 rows containing missing values (`geom_bar()`).

Residual plot:

preds<-posterior_predict(M.pe5_int.mig2,ndraws=250,summary=FALSE)

M.pe5_int.mig2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_proced$Proced_Equity_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.pe5_int.mig2.resids) #residuals look good

summary(M.pe5_int.mig2)

##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Proced_Equity_int ~ Gender * migrant2 + Education_cat3 + wealth2_sd + Age_num_sd + Marital_status + (1 | village) 
##    Data: data_fiji_proced (Number of observations: 168) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.41      0.27     0.02     1.03 1.00     2280     2620
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]         -2.97      0.49    -3.96    -2.03 1.00     3280     3171
## Intercept[2]         -1.81      0.43    -2.68    -0.97 1.00     3294     3165
## Intercept[3]         -1.25      0.41    -2.06    -0.41 1.00     3206     3249
## Intercept[4]          0.69      0.41    -0.09     1.51 1.00     3259     3175
## Genderm               0.37      0.40    -0.41     1.17 1.00     3230     2941
## migrant22            -0.45      0.39    -1.22     0.30 1.00     3496     2510
## Education_cat31      -0.16      0.41    -0.96     0.62 1.00     3043     2908
## Education_cat33      -0.88      0.37    -1.63    -0.16 1.00     3163     3091
## wealth2_sd            0.12      0.32    -0.50     0.72 1.00     2914     3030
## Age_num_sd            0.20      0.37    -0.52     0.92 1.00     3306     2933
## Marital_statuss       0.57      0.50    -0.37     1.57 1.00     3259     3256
## Marital_statusw      -0.23      0.50    -1.17     0.75 1.00     3105     2972
## Genderm:migrant22    -0.17      0.67    -1.47     1.15 1.00     3210     3285
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

loo.M.pe5_int.mig2<-loo(M.pe5_int.mig2)

Gender x marital status

#```{r M.pe5_int.marital2, cashe=TRUE} M.pe5_int.marital3<-brm(bf(Proced_Equity_int~Gender*Marital_status+Education_cat3+wealth2_sd+Age_num_sd+migrant2+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#{r} saveRDS(M.pe5_int.marital3,file="M.pe5_int.marital3.rda") #

M.pe5_int.marital3<-readRDS(file="M.pe5_int.marital3.rda")

M.pe5_int.marital2<-M.pe5_int.marital3

pars <- M.pe5_int.marital2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.pe5_int.marital2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.pe5_int.marital2)

 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)

 y_int <-  data_fiji_proced$Proced_Equity_int

ppc_bars(y_int, yrep_int)

Trace plots:

M.pe5_int.marital2%>%mcmc_trace

Chain convergence:

M.pe5_int.marital2%>% mcmc_plot(type='rhat_hist')

## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.pe5_int.marital2%>% mcmc_plot(type='neff_hist')

## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.pe5_int.marital2%>% mcmc_plot(type='acf_bar')

## Warning: Removed 84 rows containing missing values (`geom_bar()`).

Residual plot:

preds<-posterior_predict(M.pe5_int.marital2,ndraws=250,summary=FALSE)

M.pe5_int.marital2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_proced$Proced_Equity_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.pe5_int.marital2.resids) #residuals look good

summary(M.pe5_int.marital2)

##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Proced_Equity_int ~ Gender * Marital_status + Education_cat3 + wealth2_sd + Age_num_sd + migrant2 + (1 | village) 
##    Data: data_fiji_proced (Number of observations: 168) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.27     0.02     1.04 1.00     2905     2742
## 
## Population-Level Effects: 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]               -3.09      0.49    -4.10    -2.17 1.00     3157
## Intercept[2]               -1.92      0.43    -2.78    -1.10 1.00     3093
## Intercept[3]               -1.37      0.42    -2.19    -0.57 1.00     3109
## Intercept[4]                0.58      0.41    -0.20     1.38 1.00     3226
## Genderm                     0.14      0.40    -0.65     0.92 1.00     3327
## Marital_statuss             0.33      0.69    -1.00     1.68 1.00     3371
## Marital_statusw            -0.48      0.56    -1.55     0.62 1.00     3262
## Education_cat31            -0.11      0.41    -0.94     0.70 1.00     3240
## Education_cat33            -0.91      0.36    -1.64    -0.21 1.00     3110
## wealth2_sd                  0.15      0.31    -0.47     0.79 1.00     2986
## Age_num_sd                  0.22      0.36    -0.48     0.94 1.00     3421
## migrant22                  -0.54      0.34    -1.21     0.11 1.00     3176
## Genderm:Marital_statuss     0.50      0.92    -1.30     2.31 1.00     3429
## Genderm:Marital_statusw     1.08      1.04    -0.90     3.14 1.00     3298
##                         Tail_ESS
## Intercept[1]                3194
## Intercept[2]                3024
## Intercept[3]                3162
## Intercept[4]                3090
## Genderm                     3243
## Marital_statuss             3064
## Marital_statusw             3047
## Education_cat31             3169
## Education_cat33             3213
## wealth2_sd                  2806
## Age_num_sd                  3056
## migrant22                   3296
## Genderm:Marital_statuss     2980
## Genderm:Marital_statusw     3286
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

loo.M.pe5_int.marital2<-loo(M.pe5_int.marital2)

loo.M.pe5_int.marital2

## 
## Computed from 3200 by 168 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -240.5  9.2
## p_loo        17.8  1.2
## looic       480.9 18.4
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x education

#```{r M.pe5_int.edu2, cashe=TRUE} M.pe5_int.edu2<-brm(bf(Proced_Equity_int~Gender*Education_cat3+Marital_status+wealth2_sd+Age_num_sd+migrant2+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.pe5_int.edu2,file="M.pe5_int.edu2.rda")

M.pe5_int.edu2<-readRDS(file="M.pe5_int.edu2.rda")

pars <- M.pe5_int.edu2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.pe5_int.edu2%>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.pe5_int.edu2)

 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)

 y_int <- data_fiji_proced$Proced_Equity_int

ppc_bars(y_int, yrep_int)

Trace plots:

M.pe5_int.edu2%>%mcmc_trace

Chain convergence:

M.pe5_int.edu2%>% mcmc_plot(type='rhat_hist')

## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.pe5_int.edu2%>% mcmc_plot(type='neff_hist')

## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.pe5_int.edu2%>% mcmc_plot(type='acf_bar')

## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.pe5_int.edu2<- as.array(M.pe5_int.edu2)

Residual plot:

preds<-posterior_predict(M.pe5_int.edu2,ndraws=250,summary=FALSE)

M.pe5_int.edu2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_proced$Proced_Equity_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.pe5_int.edu2.resids) #residuals look good

summary(M.pe5_int.edu2)

##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Proced_Equity_int ~ Gender * Education_cat3 + Marital_status + wealth2_sd + Age_num_sd + migrant2 + (1 | village) 
##    Data: data_fiji_proced (Number of observations: 168) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.26     0.02     1.00 1.00     2776     2989
## 
## Population-Level Effects: 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]               -3.05      0.48    -4.01    -2.14 1.00     3201
## Intercept[2]               -1.89      0.43    -2.77    -1.05 1.00     3144
## Intercept[3]               -1.33      0.42    -2.16    -0.52 1.00     3236
## Intercept[4]                0.62      0.41    -0.16     1.43 1.00     3286
## Genderm                     0.24      0.44    -0.63     1.09 1.00     3284
## Education_cat31            -0.44      0.54    -1.51     0.60 1.00     3246
## Education_cat33            -0.90      0.45    -1.79    -0.02 1.00     3196
## Marital_statuss             0.53      0.50    -0.45     1.51 1.00     3352
## Marital_statusw            -0.16      0.51    -1.14     0.89 1.00     3253
## wealth2_sd                  0.08      0.32    -0.56     0.70 1.00     3169
## Age_num_sd                  0.21      0.36    -0.48     0.92 1.00     3069
## migrant22                  -0.53      0.33    -1.18     0.11 1.00     3303
## Genderm:Education_cat31     0.45      0.77    -1.06     1.96 1.00     3321
## Genderm:Education_cat33     0.01      0.70    -1.43     1.36 1.00     3329
##                         Tail_ESS
## Intercept[1]                2827
## Intercept[2]                3022
## Intercept[3]                3053
## Intercept[4]                3223
## Genderm                     3174
## Education_cat31             3031
## Education_cat33             3016
## Marital_statuss             3131
## Marital_statusw             3044
## wealth2_sd                  3171
## Age_num_sd                  3078
## migrant22                   3000
## Genderm:Education_cat31     3170
## Genderm:Education_cat33     3078
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

loo.M.pe5_int.edu2<-loo(M.pe5_int.edu2)

loo.M.pe5_int.edu2

## 
## Computed from 3200 by 168 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -241.1  9.3
## p_loo        18.4  1.2
## looic       482.2 18.7
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x age

#```{r M.pe5_int.age2, cashe=TRUE} M.pe5_int.age2<-brm(bf(Proced_Equity_int~Gender*Age_num_sd+Education_cat3+Marital_status+wealth2_sd+migrant2+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#``

#saveRDS(M.pe5_int.age2,file="M.pe5_int.age2.rda")

M.pe5_int.age2<-readRDS(file="M.pe5_int.age2.rda")

pars <- M.pe5_int.age2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.pe5_int.age2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.pe5_int.age2)

 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)

 y_int <- data_fiji_proced$Proced_Equity_int

ppc_bars(y_int, yrep_int)

Trace plots:

M.pe5_int.age2%>%mcmc_trace

Chain convergence:

M.pe5_int.age2%>% mcmc_plot(type='rhat_hist')

## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.pe5_int.age2%>% mcmc_plot(type='neff_hist')

## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.pe5_int.age2%>% mcmc_plot(type='acf_bar')

## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.pe5_int.age2<- as.array(M.pe5_int.age2)

Residual plot:

preds<-posterior_predict(M.pe5_int.age2,ndraws=250,summary=FALSE)

M.pe5_int.age2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_proced$Proced_Equity_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.pe5_int.age2.resids) #residuals look good

summary(M.pe5_int.age2)

##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Proced_Equity_int ~ Gender * Age_num_sd + Education_cat3 + Marital_status + wealth2_sd + migrant2 + (1 | village) 
##    Data: data_fiji_proced (Number of observations: 168) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.41      0.27     0.03     1.05 1.00     2690     2872
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]          -3.04      0.48    -3.99    -2.12 1.00     3087     2954
## Intercept[2]          -1.87      0.42    -2.67    -1.05 1.00     2809     2971
## Intercept[3]          -1.32      0.40    -2.09    -0.51 1.00     2981     3131
## Intercept[4]           0.63      0.40    -0.12     1.38 1.00     2962     2877
## Genderm                0.30      0.36    -0.40     1.01 1.00     3196     3054
## Age_num_sd             0.26      0.50    -0.73     1.22 1.00     2974     3053
## Education_cat31       -0.24      0.41    -1.02     0.57 1.00     3342     3290
## Education_cat33       -0.91      0.36    -1.61    -0.21 1.00     3232     3017
## Marital_statuss        0.54      0.52    -0.45     1.59 1.00     3098     2827
## Marital_statusw       -0.25      0.53    -1.31     0.82 1.00     2922     2937
## wealth2_sd             0.11      0.31    -0.54     0.69 1.00     3280     3204
## migrant22             -0.52      0.32    -1.17     0.11 1.00     3290     3089
## Genderm:Age_num_sd    -0.08      0.64    -1.31     1.18 1.00     2900     2926
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

loo.M.pe5_int.age2<-loo(M.pe5_int.age2)

loo.M.pe5_int.age2

## 
## Computed from 3200 by 168 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -240.1  9.3
## p_loo        17.1  1.1
## looic       480.2 18.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x wealth

#```{r M.de5_int.wealth2, cashe=TRUE} M.pe5_int.wealth2<-brm(bf(Proced_Equity_int~Gender*wealth2_sd+Age_num_sd+Education_cat3+Marital_status+migrant2+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.pe5_int.wealth2,file="M.pe5_int.wealth2.rda")

M.pe5_int.wealth2<-readRDS(file="M.pe5_int.wealth2.rda")

pars <- M.pe5_int.wealth2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.pe5_int.wealth2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.pe5_int.wealth2)

 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)

 y_int <- data_fiji_proced$Proced_Equity_int

ppc_bars(y_int, yrep_int)

Trace plots:

M.pe5_int.wealth2%>%mcmc_trace

Chain convergence:

M.pe5_int.wealth2%>% mcmc_plot(type='rhat_hist')

## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.pe5_int.wealth2%>% mcmc_plot(type='neff_hist')

## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.pe5_int.wealth2%>% mcmc_plot(type='acf_bar')

## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.pe5_int.wealth2<- as.array(M.pe5_int.wealth2)

Residual plot:

preds<-posterior_predict(M.pe5_int.wealth2,ndraws=250,summary=FALSE)

M.pe5_int.wealth2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_proced$Proced_Equity_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.pe5_int.wealth2.resids) #residuals look good

summary(M.pe5_int.wealth2)

##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Proced_Equity_int ~ Gender * wealth2_sd + Age_num_sd + Education_cat3 + Marital_status + migrant2 + (1 | village) 
##    Data: data_fiji_proced (Number of observations: 168) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.27     0.02     1.03 1.00     2777     2744
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]          -3.04      0.47    -3.97    -2.12 1.00     2865     2963
## Intercept[2]          -1.87      0.42    -2.70    -1.05 1.00     3069     3021
## Intercept[3]          -1.31      0.40    -2.10    -0.51 1.00     3137     2842
## Intercept[4]           0.63      0.40    -0.13     1.46 1.00     3284     3150
## Genderm                0.31      0.37    -0.41     1.04 1.00     3132     2800
## wealth2_sd             0.16      0.39    -0.61     0.94 1.00     2820     2527
## Age_num_sd             0.22      0.37    -0.51     0.95 1.00     3135     3091
## Education_cat31       -0.24      0.40    -1.03     0.57 1.00     3047     2839
## Education_cat33       -0.91      0.37    -1.65    -0.21 1.00     3074     2985
## Marital_statuss        0.55      0.52    -0.46     1.55 1.00     3234     2894
## Marital_statusw       -0.23      0.49    -1.21     0.70 1.00     3185     3064
## migrant22             -0.52      0.33    -1.16     0.12 1.00     3134     3241
## Genderm:wealth2_sd    -0.18      0.62    -1.40     1.03 1.00     3144     3004
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

loo.M.pe5_int.wealth2<-loo(M.pe5_int.wealth2)

loo.M.pe5_int.wealth2

## 
## Computed from 3200 by 168 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -240.1  9.3
## p_loo        17.0  1.1
## looic       480.1 18.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Compare models

loo_compare(loo.M.pe5_int.wealth2,loo.M.pe5_int.age2,loo.M.pe5_int.edu2,loo.M.pe5_int.marital2,loo.M.pe5_int.mig2,loo.M.pe5_add_2 )

##                    elpd_diff se_diff
## M.pe5_add_2         0.0       0.0   
## M.pe5_int.wealth2  -1.0       0.3   
## M.pe5_int.age2     -1.1       0.2   
## M.pe5_int.mig2     -1.4       0.4   
## M.pe5_int.marital2 -1.4       1.1   
## M.pe5_int.edu2     -2.0       0.7